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ABSTRACT 

The finite- volume and finite-difference implementations of high-order accurate essentially 
non-oscillatory shock-capturing schemes are discussed and compared. Results obtained with 
fourth-order accurate algorithms based on both formulations are examined for accuracy, 
sensitivity to grid irregularities, resolution of waves that are oblique to the mesh, and com- 
putational efficiency. Some algorithm modifications that may be required for a given appli- 
cation are suggested. Conclusions that pertain to the relative merits of both formulations 
are drawn, and some circumstances for which each might be useful are noted. 
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Introduction 

Two basic forimilations, finite-voliime and finite-difference, for the implemen- 
tation of high-order accurate, essentially non-oscillatory (ENO) shock-capturing 
schemes have been the subject of considerable interest in recent years. These 
schemes achieve high-order spatial accuracy in smooth regions by means of a piece- 
wise polynomial approximation operator that is also designed to avoid oscillations 
associated with interpolation across steep gradients. As such, they are well suited 
for the study of aeroacoustic and transition-related problems and may serve as an 
alternative to spectral methods for solving such problems when shocks or complex 
geometries are involved. 

The finite-voliime implementation, first presented by Harten et al.^, is preferred 
for its strict adherence to the integral form in which conservation laws are defined. 
The primary motivation for the use of the finite-difference approach of Shu and 
Osher^ is computational efficiency. These formulations are briefly described, after 
which results of their numerical implementations are presented for comparison. 
The intent of this work is to acquaint the reader with the relative merits of both 
formulations, the circumstances for which each might be useful, and some details of 
implementation that may be required for a given application. The performances 
of both algorithms are compared for accuracy, sensitivity to grid irregularities, 
resolution of waves that are oblique to the mesh, and computational efficiency. 

Discrete Formulations 

The finite-volume and finite-difference algorithms compared in this paper differ 
fimdamentally in the way a system of equations is solved. In both cases, a weak 
solution of a system of conservation laws is ultimately obtained. The conservation 
of some quantity U in a. spatial domain D can be written 

/ UdV = - f F ‘H dS (1) 

dt Jd JdD 

where F is the flux, dD is the boundary of D, dV is a volume element of D, dS is an 
element of surface area on dD, and n denotes the outward imit normal to dD. 

In the finite-volume approach, the conservation law itself in Eq. (1) is approxi- 
mated. The spatial domain is discretized, D = {Di}, which results in 

= -V f 

Ot Vi Jojj. 
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where Vi is the volume of Di and 

tJi = ^ f U dV (3) 

JDi 

is the cell average of U in D,. Eq. (2) is solved for all i, which yields a solution {[/, } 
of cell averages. 

In the finite-difference approach, a pointwise solution is desired. To this end, time 
differentiation and spatial integration are interchanged in Eq. (2), the divergence 
theorem is applied on the right-hand side and, in the limit as Vi —>■ 0, 

= -(vf). (4) 

Special care must be taken when this formulation is implemented, because flux 
conservation is not as readily achieved as in Eq. (2). 

For integration in time, the method of lines will be employed for Eqs. (2) and 
(4). High-order accurate Runge-Kutta methods, developed by Shu and Osher^, are 
implemented in the finite-volume and finite-difference schemes to be compared in 
this paper. Hence, the brief description of the schemes to follow will concern only 
the right-hand sides of Eqs. (2) and (4). 

Both discrete algorithms involve a reconstruction step followed by an evolution 
step. What is meant by reconstruction is a high-order accurate polynomial approx- 
imation at some point in time. In the finite-volume formulation, the solution U is 
reconstructed from the cell averages to high order within each cell Di and evaluated 
on the boundary 5Z>, . The evolution step involves the solutions of the local Riemann 
problems that arise from the piecewise continuous reconstruction (See, e.g.. Refs. 1 
and 3 and the references therein.). The spatial integration along the boimdary dDi is 
achieved by a correspondingly high-order quadrature. This method will be referred 
to throughout the paper as ENO-FV. Shu and Osher^>^ have proposed the use of 
the finite-difference formulation in Eq. (4), in which the reconstruction operator is 
applied directly to the pointwise flux values. The evolution step arises from a flux- 
vector splitting strategy that is built into the reconstruction. The term ENO-FD 
will be used to refer to this finite-difference scheme. The pointwise nature of this 
formulation eliminates the need for dealing with cell averages or the integration of 
a flux over the boimdary of a cell. These distinctions between the two formulations 
become most important with regard to the issue of cost, and will be discussed in 
more detail in a later section. 
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The most unique aspect of ENO reconstruction operators is their use of adaptive 
stenciling. At the i-th discrete location, the reconstruction set that determines a 
pol 3 momial Vi is chosen by a searching algorithm in which decisions are based 
on local smoothness criteria. For example, in one dimension, where the {D,} are 
intervals on the Real line, Harten et al.^ have suggested a hierarchical algorithm 
that is based on local sets of A;-th divided differences where i - m < j < i + m, 

A: = 1,2, . . . ,m, and m is the desired degree of Vi. If a contiguous one-dimensional 
stencil for Vi is defined by its left-most index im, then the first-order choice is clearly 
ii = i. Then, for fc = 2, 

i = / ’ if \^ik-i\ < \^i\ (5) 

\ ik-i , otherwise 

Because this algorithm allows the reconstruction stencil to shift freely with the 
detection of any numerical gradient, it will be referred to as *Treely adaptive. 
Reasons for modif 5 dng this search will become apparent. 

A few observations that concern the available t 3 rpes of reconstruction operators 
are in order. Within the Shu-Osher approach, when conservation is desired, the 
implementation of the high-order reconstruction operator requires a uniform compu- 
tational mesh. Therefore, the application of the ENO-FD algorithm on a non-uniform 
physical domain requires a sufficiently smooth transformation to a uniform mesh if 
third-order or higher accuracy is desired. An analogous form of this transformational 
reconstruction (TR) can also be implemented within the finite-volume formulation. In 
this case, the reconstruction operator is applied to a set of volume-weighted averages 
{ViUi}, and therefore sufficient smoothness in a spatial transformation is required.® 
However, another option exists for the ENO-FV algorithm. The polynomial approxi- 
mation can be performed in physical space, which releases such burdensome restric- 
tions on grid smoothness.®’^’® This latter procedure will be referred to as physical 
reconstruction (PR). Although this PR operator poses no problems in one dimension, 
its implementation in multidimensional space can be quite complex when the imple- 
mentation must allow for local stencil adaptation,'^’® However, the multidimensional, 
finite-volume TR operator can be readily implemented because it is defined as a prod- 
uct of one-dimensional operators.® When required, the finite-volume algorithms will 
be distinguished as ENO-FV-TR and ENO-FV-PR. 

One-Dimensional Rarefaction Wave 

The first test case involves the solution of the Euler equations of gas dynamics in 
one spatial dimension, as it pertains to the movement of a right-traveling rarefaction 
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wave through a domain of highly varying mesh spacing. The effects of two mesh 
transformations will be examined. One transformation has the required smoothness 
for a fourth-order TR operator and the other does not. 

For both transformations, the uniform computational domain is given by 
{ ^ i ^ 6 } . The subset { —4 < ^ < 4 } is divided into five intervals: 
[—4, —3], [—3, —1], [—1, 1], [1,3], and [3,4] . In each of these interior intervals, the map- 
ping X = x((f) causes the physical mesh width Ax to vary rapidly; the intervals 
{ — 6 < ^ < -4 } and { 4 < ^ < 6 } are mapped uniformly. 

In the fc-th interior interval, the first transformation is of the form 


a:(0 = 6 + sin 




(6) 


where is an element of the set {-4, -2, 0, 2, 4} . The mesh spacing Ax in the uniform 
regions (fore and aft) is determined such that the connections at ^ = ±4 are smooth. 

To generate a smoother gnd on { —4 < ^ ^ 4 } , a mapping of the following form 
is used: 


x{i) = 






15 3 1 

— sin (tt^) sin (27 t^) + — sin 

7T 7T 3tT 


(7) 


The parameters a and ^ are determined so that the ratio of the maximum to 
minimum values of x^ are identical for the two grids, and that the physical distance 
between x(— 4) and x(4) is the same for both. 

Fig. 1(a) illustrates the similar behavior of these two transformations on 
0 < ^ < 2 . The seemingly odd formula in Eq. (7) was chosen because its derivative 
is of the form 


x^ = a-^ sin® (| ^ ) 

The transformation derivatives are plotted on the same interval in Fig. 1(b). The 
value .f^ represents the derivative normalized by its maximum value. For the 
transformation (6), x^^ is discontinuous for ^ = dbl,±3, which makes for an overall 
gnd. The mesh produced by Eq. (7) will be C°° on the nonimiform region and C® 
overall, because of the connections with the imiform intervals at ^ = ±4 . 

The initial solution consists of an isentropic expansion, smoothly distributed on 
—6.5 < X < —5.5, with constant states, U\ and U-i, on the right and left, respectively. 
The strength of this rarefaction is determined by requiring a mean temperature 
change of ±5 percent across the wave and supersonic flow on either side. The problem 
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is nondimensionalized with respect to the mean solution. For f > 0, the rarefaction 
wave moves to the right, and its behavior is monitored within the nonuniform region 
of the mesh until t = 4.0. 

Both algorithms used here are fourth-order accurate in the L\ sense. Fig. 2 
compares solutions for the two formulations by x-t contour plots of the density, for 
0.0 < < < 4.0 , on a mesh of 120 intervals with the transformation given by Eq. (6). 
The freely adaptive stencil algorithm in Eq. (5) is employed. All of the contour plots 
in this section represent the entire solution range, pi < p < P 2 , with 20 equally spaced 
levels. The ENO-FD solution exhibits significant qualitative error in neighborhoods 
about the second-derivative discontinuities in the mesh. In Fig. 2(b), such problems 
are not apparent when the solution is computed with the freely adaptive ENO-FV- 
PR algorithm. 

Fig. 3 depicts the use of the same fourth-order ENO-FD algorithm on the 
grid generated by Eq. (7). With the required mesh smoothness to support the 
reconstruction, no visible distortions appear as on the grid. However, the fourth- 
order design accuracy is not achieved for this solution when density error is measured 
with respect to the Li norm in a region to the left of the rarefaction att = 4.0 . The 
same accuracy problem was found with the freely adaptive ENO-FV-PR scheme. The 
results of both mesh-refinement studies are shown on a log-log plot in Fig. 4(a). The 
spatial discretizations employed in this study are 120, 240, 480, and 960 intervals. 
The numbers on the loci represent the computational order of accuracy as measured 
between the two finest meshes. 

A similar loss-of-accuracy phenomenon was reported by Rogerson and Meiberg®, 
which prompted a response from Shu^® that the problem arises by allowing the 
stencil to adapt too freely. Shu has suggested that the stencil adaptation algorithm 
be modified to bias the stencil towards one that is stable, in the sense of linear 
stability analysis. In the present application, the resulting stencil is one which is 
upwind biased. This biasing can be done by implementing a factor <r in the stencil 
search in Eq. (5), vis. 

zjt = I ^ I < (9) 

\ U_i , otherwise 

where {oi^ctr) = (1,5-) or (a, 1), for biasing to the left or right, respectively, with 
a > \, 

Fig. 4(b) shows grid-refinement results that are analogous to those in Fig. 4(a) 
using the modification in Eq. (8) with a = 2.0; however, a fourth-order error reduction 
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is still not evident because the error is measured in a region where numerical 
gradient? are extremely small and ratios of neighboring gradients may be much 
larger than the chosen biasing parameter. Atkins^^ has suggested that an additional 
constraint on the stencil adaptation is required. The purpose of this constraint is 
to bias the stencil toward one that is stable wherever the solution is smooth. A 
parameter installed for this purpose can be considered a lower threshold for the 
magnitudes of the local differences {6j^} , below which a stencil is forced toward 
a stable target, regardless of the relative magnitudes of neighboring numerical 
gradients. The present implementation of such a parameter can be written 


if 



1 


< e 


and 



< e 


then 


h- i*k ) 


( 10 ) 


where e is a small parameter, and defines the stable stencil at the k-th level. 

Another grid-refinement study was performed, with Eqs. (8) and (9) wherein 
a = 2.0 and e = 0.01. The results are shown in Fig. 5(a) and are much more consistent 
with the scheme’s design accuracy. In fact, the ENO-FV-PR algorithm, coupled with 
the stenciling modifications in Eqs. (8) and (9), 5 delds fourth-order computational ac- 
curacy on the grid, which is shown in Fig. 5(b). The second-order convergence 
exhibited by the ENO-FD algorithm in this plot is expected. However, Fig. 6 shows 
that the qualitative error in the ENO-FD solution on the grid, that was previ- 
ously shown in Fig. 2(a), has been significantly reduced with these stencil biasing 
modifications. 

At this point, the apparently greater robustness of the ENO-FV algorithm can be 
attributed to the fact that this formulation allows for a physical reconstruction that 
is not available for the finite-difference algorithm. However, as previously noted, 
multidimensional extensions for this more generalized operator can be complicated 
and costly. Therefore, some results produced with the ENO-FV-TR algorithm will be 
discussed. In Fig. 7(a), when the freely adaptive version of this algorithm is imple- 
mented on the C® grid, significant error exists that was not evident in the analogous 
case for the ENO-FD solution in Fig. 3. This error is due to the reconstruction of 
volume-weighted averages, in which a rapidly varying mesh will have an inordinate 
effect on the stencil choice in a region of small solution gradients. However, biasing 
the reconstruction stencils with Eqs. (8) and (9) is not enough to entirely rid the 
solution of these grid distortion errors, as shown in Fig. 7(b). The problem is still 
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related to the reconstruction of volume-weighted averages, in particular to the fact 
that a imiform flow is not preserved by this procedure if the anal 5 rtic transformation 
is used for the necessary discrete values . Therefore, instead of using Eq. (7b) to 
compute the derivatives, these derivatives are numerically approximated in a man- 
ner that enables the reconstruction operator to preserve the free stream exactly®. 
The resulting solution, with this final modification, is shown in Fig. 7(c). In this 
final form, the numerical accuracy of the ENO-FV-TR algorithm is foimd to perform 
to design as shown in Fig. 8(a). The ENO-FD and ENO-FV-PR results from Fig. 5(a) 
are repeated in Fig. 8(b), for comparison. 

Two-Dimensional Channel Flow 

The ENO-FD and ENO-FV algorithms are now compared in two spatial dimen- 
sions. The test case involves a steady, subsonic flow in a channel of var 3 dng area. 
Although high-order ENO schemes are clearly designed with unsteady solutions in 
mind, many such solutions of interest can be considered as the imposition of per- 
turbations upon a steady flow. It is therefore essential that a steady flow be accu- 
rately predicted in order to obtain meaningful results from imsteady problems of an 
aeroacoustic or transitional nature. Unless otherwise stated, the remainder of the 
applications of the ENO-FV algorithm will involve the TR type of reconstruction. 

This channel flow solution is assumed to be governed by the two-dimensional 
Euler equations and is computed with both the ENO-FD and ENO-FV algorithms 
on two different geometries. For the two channels imder consideration, the length- 
to-height ratio is LjH = 1.5, with constant-area sections fore and aft of a section 
of varying area. Each constant-area section has length L/5; at the throat, the con- 
striction is 10 percent of H. The most significant difference between the two geome- 
tries is the order of smoothness to which the constant-area sections are connected 
to the center section. Both channels can be described as follows. Let the rectangle 
{0 < ^ < L} X {0 < t] < H} denote the computational domain for each channel. The 
identity maps (^,r/) to the physical point (x,y ) , for ^ in [0, L/5] U [4L/5,L] and all r/. 
On the interval {L/5 < ^ < 4L/5} , the geometry for the varying-area section is given 
by a transformation of the form 

X = ^ (11a) 

y = ]jyz{() + {i-jj)vm (lib) 

where yi{i) and j/2(0 respectively, the equations for the bottom and top walls. 
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The first middle section considered is determined by walls of the form 

2/i(0 = a« sin^ (6,-^ + c.) + </,■, i = 1,2 (12) 

where the values of the coefficients are set for the desired throat constriction and 
a smooth connection to the constant-area sections. Three continuous ( derivatives 
of y exist at ^ = Lfb^ALfb, for all t/. Therefore, Eqs. (10) and (11) generate a 
geometry. Fig. 9 illustrates this geometry on a 60x40 mesh. The second channel 
geometry differs from the first only in that the walls in the middle section are given 
by polynomials of the form 

• yiiO = + Ci^\ i = 1,2 (13) 

In this case, the connections of the middle to the outer sections are continuous in ^ 
to only one derivative, for 0 < y < H . 

The desired solution is that of a steady-state flow that is caused by a imiform, 
parallel free stream entering the channel at a; = 0 . This solution is achieved by 
solving the time-dependent Euler equations with the fourth-order algorithms. At 
^ = 0 , the solution in the entire channel is set to free stream with a Mach number of 
0.3. Tangency is imposed on the walls, and a non-reflecting boundary condition^^ is 
applied at the inflow and outflow. The stencil-biasing modifications in Eqs. (8) and (9) 
are essential for convergence of the solution to a steady state, and are implemented 
with a = 2.0 and e = 0.01 . Flux residuals were readily driven to machine zero in 
both test cases. 

These solutions were computed by both algorithms on a sequence of successively 
refined grids, and the solutions’ errors were determined by deviation from isentropy, 
as measured by the quantity S = P/p'^ . The refinement sequence employed here 
is 15x10, 30x20, 45x30, 60x40, and 90x60. The results of the grid-refinement 
study on the geometry are shown in Figs. 10(a), (b),and (c). The “Global” er- 
ror is computed over the entire computational domain, the “Wall” error is com- 
puted only at the points along one wall, and the “Interior” error is computed on 
{L/A < ( < SL/A} U {H/A <T} < 3/f/4) . Both algorithms evidently perform at or near 
design accuracy, as expected for fourth-order schemes on a mesh. 

In Figs. 10(d), (e), and (f) are the grid refinement results for the mesh. 
As expected, second-order results are obtained on the wall for both algorithms. 
However, the finite-volume algorithm performs at third-order accuracy with respect 
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to the global error and at design accuracy on the interior, while the finite-difference 
algorithm shows second-order accuracy in all three measures. These results suggest 
that, for the ENO-FD algorithm, the second-order entropy error that arises from 
the non-smoothness at the section connections is propagating into the interior. This 
assumption is supported by Fig. 11, in which the quantity log 5 is plotted along the 
center line of a 60x40 mesh. Clearly, in the finite-difference solution, there is a 
lower-order entropy error within the middle section than exists in the constant-area 
sections. 

Oblique Sod’s Problem 

The final test imder consideration is selected to compare the capabilities of the 
two algorithms to resolve waves that are oblique to the computational mesh. Sod’s 
problem^^ will be solved in two-dimensional space so that the planar waves produced 
will propagate at various angles of incidence with respect to a rectangular gnd. The 
intent is not only to inspect the qualitative resolution of the oblique waves, but also 
to quantify the manner in which each algorithm detects an oblique wave with respect 
to its detection of a wave that is normal to the mesh. 

Sod’s problem is a Riemann problem which is subject to the solution of the Euler 
equations and determined by initial conditions that consist of a thermodynamic 
discontinuity imposed upon a fluid at rest. The magmtudes of the density and 
pressure jumps are, from left to right, plIpr = 8 and Pl/Pr = 10, respectively. 
Such an initial solution is installed on a two-dimensional Cartesian grid so that the 
discontinuity is a straight line that makes an angle 6 with the x axis. The physical 
mesh varies with 0 in the following manner. A computational mesh is defined on 
the rectangle [0,1] x [0,i/j, with a length-to-height ratio LfH = &. If 0 < 0 < 1^12, 
then the physical domain is {0 < x < L 0 ] x {0 < y < Hq} , and is determined by 
He = HlsinO and Le = I/sin^. This scaling achieves the same grid resolution 
normal to the wave propagation on a gdven mesh at some fixed time for all choices 
of 0. The mesh discretization in each direction is chosen such that Ar = Ay . 

At i = 0, the initial discontinuity is positioned at (a:,y) = (3L/8,0) and inclined 
at the angle 6. The angles of inclination are chosen so that tan(? is an integer and 
that the boundary conditions at y = 0 and y = H$ can be determined in a “shifted- 
periodic” manner. In particular, the test angles are 0 = arctan 1, arctan 2, and arctan 4 , 
in addition to the one-dimensional problem 0 = 1:12 . The solution is computed to 
t = \. 2 The stencil biasing in Eq. (8) is used in the reconstruction procedures with 
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a = 2.0. 

Fig. 12(a) represents the solution at < = 1.2, on a 96x16 grid, with 0 = arctan 1, 
with the use of the ENO-FV algorithm. The three wave structures are, from left to 
right, an expansion wave, a contact discontinuity, and a shock. The axis variables x' 
and y' represent the physical coordinates scaled by sin 6 . Fig. 12(b) depicts the same 
solution as in Fig. 12(a), for both algorithms, along the grid line y' = 0. On this level, 
the qualitative difference between the solutions is barely detectable. Therefore, in 
Fig. 13, the solutions are examined in the vicinity of the contact discontinuity for 
all test values of 6. (The successive loci in this plot have been shifted upward and 
to the right.) Although the ENO-FV scheme produces a slightly steeper numerical 
gradient at the location of the discontinuity, such a small difference can be judged 
insignificant. Fig. 13, then, suggests that both algorithms perform equally well 
with respect to oblique waves irrespective of 9. On a more quantitative level, 
the differences between the oblique cases and the one-dimensional case, for each 
algorithm, have been measured. These results are plotted in Figs. 14(a) and 14(b). 
By this measure, the finite-volume scheme performs marginally better, particularly 
with respect to the rarefaction wave. 

Cost Comparison 

As with any numerical algorithm, the cost of implementing either the ENO-FD 
or ENO-FV scheme is a major concern. Conceptually, the two formulations can be 
made equally cost effective for high-order solutions of one dimensional problems or for 
first or second-order accurate solutions in multiple dimensions. However, for third 
or higher-order accuracy in two or more dimensions, the algorithms are radically 
different. This difference translates to a significant disparity in cost. 

With regard to cost, the ENO-FD algorithm has a clear advantage when applied 
to multidimensional problems. This advantage can be entirely attributed to the fact 
that the finite-difference operator solves a system of equations in a pointwise man- 
ner. In this case, a high-order multidimensional reconstruction can be accomplished 
in a dimension-by-dimension fashion. If the algorithm is r-th order accurate in the 
L\ sense, then the reconstruction stencil in k dimensions contains k{r - 1) -f 1 points. 
The dimension-by-dimension approach to reconstruction is not possible within the 
finite-volume formulation, if third-order or better accuracy is required. Because the 
solution at any time is in the form of cell averages, the multidimensional ENO-FV- 
TR reconstruction is implemented as a product of one-dimensional operators, and 
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thereby requires cells in a reconstruction stencil. Moreover, the Gaussian quadra- 
ture used to integrate the flux on each cell boundary requires multiple solutions 
of the Riemann problem on each cell interface. The number of points n, required 
in this quadrature on an interface is Uq = [int((r+l)/2)]*^“^ , where int(-) denotes 
integer truncation. This high-order quadrature is entirely avoided within the finite- 
difference approach. It should be noted that both formulations require the same 
amoimt of logic in the stencil adaptation, which is the most expensive part of either 
algorithm. 

Fig. 15 depicts the relative CPU cost of the two algorithms, for solutions of 
the Euler equations in two dimensions, for 2nd, 3rd, 4th, and 5th-order accurate 
applications. CPU times were measured on a Cray-YMP and normalized with 
respect to the 2nd-order ENO-FD cost. An estimated cost comparison for three- 
dimensional problems can be obtained by the following reasoning. For a specified 
order of accuracy, let Ci denote the cost of solving one equation at a point in one 
dimension; this is the same for either algorithm. Then, because of its dimension-by- 
dimension implementation, the estimated cost incurred by the ENO-FD algorithm 
in solving the Euler equations in k dimensions on a grid of points is 

CFD = Ci{k + 2)kN^ (14) 

This formula predicts a cost-per-point increase factor of 8/3 when the one-dimensional 
ENO-FD algorithm is extended to two dimensions. This value has been supported 
by computer measurements. The corresponding factor for the extension from two to 
three dimensions is 15/8. The estimation of a three-dimensional cost for the ENO- 
FV algorithm cannot be done by such a simple linear extrapolation. There are two 
significant costs in addition to the base value given by Eq. (13). One is the extension 
of the finite-volume reconstruction to a three-dimensional product. The other arises 
from the additional 3 (n, - 1) flux computations that are required in the integration 
over the surface of each cell. These additional costs were extrapolated from two- 
dimensional CPU measurements on a Cray-YMP. Fig. 16 represents the estimated 
three-dimensional cost comparison where, again, the data are normalized by the 
2nd-order cost. 

Concluding Remarks 

For accuracy on a sufficiently smooth mesh and resolution of oblique waves, both 
algorithms perform equally well. The finite-volume implementation is less sensitive 
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to derivative discontinuities, whether in the computational mesh or in the solution. 
In particular, the ENO-FV-PR algorithm has the capacity to perform at design 
accuracy, independent of the mesh. Although the generalized multidimensional 
adaptive-stenciling implementation of this algorithm has not sufficiently matured, 
some promising results can be found in Ref. 8. Either of the multidimensional finite- 
volume algorithms is significantly more costly than the finite-difference algorithm. 
Therefore, if, for a given application, the computational domain is known to be 
sufficiently smooth and can be suitably structured, then the ENO-FD algorithm is 
the method of choice. However, for problems with complex geometries, it might pay 
to use the more expensive algorithm if the grid is significantly less costly to generate 
in a less restrictive fashion. 
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Fig. 7. Density, on the C® mesh, ENO-FV-TR. 
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Fig. 8. Li density error, t = 4.0, grid, (a) Biased stencil, free-stream preserving, (b) Biased stencil 
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Fig. 10. Li entropy errors, steady 2D channel flow. 
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Fig. 14. (a) Density deviation front ID, ENO-FD. (b) Density deviation from ID, ENO-FV. 
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